#-------------------------------------------------------------------------------
# Name:       programme batch
# Purpose:
#
# Author:      P. Sanchez, M. Carrere
#
# Created:     07/07/2011
# Copyright:   (c) BDI 2011
# Licence:     GNU
#-------------------------------------------------------------------------------
#!/usr/bin/env python
import numpy as np
import pylab as plt
import scipy as sc
from numpy import sin
from scipy import pi
import os.path
from optparse import OptionParser

from modreststate import *

import os
import sys
import guidata

def clearall():
    all = [var for var in globals() if var[0] != "_"]
    for var in all:
        del globals()[var]


parser = OptionParser()
(options, args) = parser.parse_args()
fname=str(args[0])
app = guidata.qapplication()
plt.figure(1)

dtbeforeTMS=10*1e-3
deltaBeforeTMS=500*1e-3
deltaAfterTMS=1200*1e-3
dtAfterTMS=10*1e-3
nbptFFT=2048



path="C:\ProgsPython\pyreststate\mat"  # insert the path to the directory of interest here

print fname
filename = path  + "\\" + fname
nouveau_nom = '.'.join(fname.split('.')[:-1])
print "nouveau nom = ", nouveau_nom
#exit(0)
voies=[]
del(voies)
nbchannel,nbblock,blocksize,voies,srate, dt=mat2pydata(filename)
print " 1 bnb channel=", nbchannel,"nb block=", nbblock, "bock size=", blocksize
#-----------------------------------------------------------------------------
#   axe  des temps TMS au centre
#-----------------------------------------------------------------------------
temps=[]
t_sav=[]
temps=CreatTimeTMS(blocksize,srate)
t_sav=temps                          # sauvegarde axe des temps

#-----------------------------------------------------------------------------
#   selection des voies
#-----------------------------------------------------------------------------
for numv in range(1,nbchannel+1):
    print "numero voie = ",numv
    numvoie=numv-1
    v=voies[numvoie,:,:]  # recupere la voie Fz
    v=v.transpose()
#-----------------------------------------------------------------------------
#     velection block avant TMS
#-----------------------------------------------------------------------------
    M_b=[]
    M_a=[]
    t_a=[]
    t_b=[]
    M_b,t_b, ind0_b, ind1_b, t0_b, t1_b=BlockBeforeTMS(v,blocksize,srate,temps,dtbeforeTMS,deltaBeforeTMS)
    M_a,t_a, ind0_a, ind1_a, t0_a, t1_a=BlockAfterTMS(v,blocksize,srate,temps,dtAfterTMS,deltaAfterTMS)
#-----------------------------------------------------------------------------
#     calcul des coefficients du filtre
#-----------------------------------------------------------------------------
    AlphaMin=8
    AlphaMax=14

    bn=coefPassBand(AlphaMin, AlphaMax, srate)


#-----------------------------------------------------------------------------
#            initialisations ondelette
#-----------------------------------------------------------------------------
    f_start=5.
    f_stop=50.
    deltafreq = 0.5
    sampling_rate = 2500.
    sizet_b=len(t_b)
    sizet_a=len(t_a)
    duree_b=sizet_b/sampling_rate
    duree_a=sizet_a/sampling_rate
    f0=2.
    scales = f0/np.arange(f_start,f_stop,deltafreq)*sampling_rate
    U=len(scales)
    V_b= sizet_b
    V_a= sizet_a
    WSUM_b=[]
    WSUM_a=[]
    WSUM_b=np.zeros((U, V_b), dtype=np.float)
    WSUM_a=np.zeros((U, V_a), dtype=np.float)
#-----------------------------------------------------------------------------
#  debut boucle sur ensemble des blocs
#-----------------------------------------------------------------------------
    for loop in range(nbblock):
        print "numero block = ",loop
#-----------------------------------------------------------------------------
#   un block de la vopython batchreststate.py ie selectionn?e
#-----------------------------------------------------------------------------
        signal_b=M_b[:,loop]
        signal_a=M_a[:,loop]
#-----------------------------------------------------------------------------
#   selection filtre
#-----------------------------------------------------------------------------

        filtered_b = []
        filtered_b.append(sc.convolve(bn, signal_b))	# filter the signal by convolving the signal with filter coefficients
        for i in xrange(0, len(filtered_b)):
               y_p_b = filtered_b[i]
        y_p_b=np.array(y_p_b)

        sig_b=y_p_b[0:sizet_b]

        filtered_a = []
        filtered_a.append(sc.convolve(bn, signal_a))	# filter the signal by convolving the signal with filter coefficients
        for i in xrange(0, len(filtered_a)):
               y_p_a = filtered_a[i]
        y_p_a=np.array(y_p_a)

        sig_a=y_p_a[0:sizet_a]
#-----------------------------------------------------------------------------
#   ondelette de Morlet
#-----------------------------------------------------------------------------
        wt_b=[]
        wt_b=compute_morlet_scalogram(sig_b,
                    f_start=f_start,
                    f_stop=f_stop,
                    deltafreq = deltafreq,
                    sampling_rate = sampling_rate,
                    f0=2,
                    normalisation = 0.,
                    wf=None
                    )
        WM_b=abs(wt_b) .transpose()
        WSUM_b=WSUM_b+WM_b
        wt_a=[]
        wt_a=compute_morlet_scalogram(sig_a,
                    f_start=f_start,
                    f_stop=f_stop,
                    deltafreq = deltafreq,
                    sampling_rate = sampling_rate,
                    f0=2,
                    normalisation = 0.,
                    wf=None
                    )
        WM_a=abs(wt_a) .transpose()
        WSUM_a=WSUM_a+WM_a
#-----------------------------------------------------------------------------
#   calcul de la moyenne des blocs
#-----------------------------------------------------------------------------
    WSUM_b=WSUM_b/(nbblock+1.)
    WSUM_a=WSUM_a/(nbblock+1.)
#-----------------------------------------------------------------------------
#  creation palette
#-----------------------------------------------------------------------------

#-----------------------------------------------------------------------------
#   affiche moyenne des analyses en ondelette
#-----------------------------------------------------------------------------
    ax1=plt.subplot(121)
    imb = plt.imshow(WSUM_b,
                                    interpolation='nearest',
                                    extent=(t0_b, t1_b, f_start, f_stop),
                                    origin ='lower' ,
                                    aspect = 'normal')
    for tick in ax1.xaxis.get_major_ticks():
                tick.label.set_fontsize(9)
    for tick in ax1.yaxis.get_major_ticks():
                tick.label.set_fontsize(9)
    imageName=nouveau_nom + "_" + str(numv)
    plt.title(imageName,fontsize=9, color='k')
    plt.xlabel('time [s]',fontsize=9, color='k')
    plt.ylabel('Freq. [Hz]',fontsize=9, color='k')
    ax2=plt.subplot(122, sharey=ax1)
    ima = plt.imshow(WSUM_a,
                                    interpolation='nearest',
                                    extent=(t0_a, t1_a, f_start, f_stop),
                                    origin ='lower' ,
                                    aspect = 'normal')
    for tick in ax2.xaxis.get_major_ticks():
                tick.label.set_fontsize(9)
    for tick in ax2.yaxis.get_major_ticks():
                tick.label.set_fontsize(9)
    plt.title("WT after TMS",fontsize=9, color='k')
    plt.xlabel('time [s]',fontsize=9, color='k')
    #plt.ylabel('Freq. [Hz]')
    #plt.setp(ax2.yticklabels, visible=False)
    plt.subplots_adjust(hspace=0.05)

    plt.show()
    plt.draw()

    plt.savefig(imageName)
    print "save ",imageName
    plt.close(1)
app.exec_()


